/*--- sqrtmain.c -------------------------- Listing 8-21 -------
 * Uses the longmath routines to perform Newton's algorithm for
 * finding square roots. Compares the result to the C floating-
 * point function sqrt(). This implementation keeps approximating
 * until Newton's method no longer generates a difference
 * or until 50 iterations have been performed. Beyond 50
 * iterations, it is assumed that the algorithm is ping-ponging
 * around the final number. Square root of 525 with MAX_SIZE
 * = 20 is an example of this problem.
 *-------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "longmath.h"

main(int argc, char *argv[])
{
    int i;
    char buffer[2 * MAX_SIZE];

    struct TermData *square;    /* the original number */

    struct TermData *xn,        /* variables in equation */
                    *xn1,
                    *product;

    struct TermData *two,       /* constants */
                    *three;

    struct TermData *temp1,     /* intermediate results */
                    *temp2,
                    *temp3,
                    *z;

    float fsquare, fsqrt;       /* the floats we'll compare to */

    square = TermCreate();
    xn = TermCreate();
    xn1 = TermCreate();
    product = TermCreate();

    two = TermCreate();
    three = TermCreate();

    temp1 = TermCreate();
    temp2 = TermCreate();
    temp3 = TermCreate();
    z = TermCreate();


    if (argc < 2) {
        printf("Usage: sqrtmain number\n");
        return;
    }

    AsciiToTerm(argv[1], square);

    /* no imaginary numbers in this version */
    if (square->sign == -1) {
        printf("Square cannot be negative\n");
        return;
    }

    /* square root of zero is zero */
    if (square->places_before + square->places_after == 0) {
        printf("0\n");
        return;
    }

    /*
     * We start with a guess of dividing square by 3, hence
     *           xn = square / 3.0
     */
    AsciiToTerm("3.0", three);

    if (!ComputeResult(square, '/', three, xn)) {
        printf("error 1\n");
        return;
    }


    printf("Square root of %s\n", argv[1]);
    TermToAscii(xn, buffer, NORMAL);
    printf("First guess is %s\n", buffer);

    AsciiToTerm("2.0", two);

    for (i = 0; ; i++) { /* repeat until root or max 50 loops */
        /* Newton's approach keeps approximating using
         * this formula:
         *  xn1 = (square + (xn * xn)) / (2.0 * xn)
         */
        TermCopy(z, xn);
        if (!ComputeResult(z, '*', xn, temp2)) {
            printf("error 2\n");
            return;
        }

        if (!ComputeResult(square, '+', temp2, temp1)) {
            printf("error 3\n");
            return;
        }

        if (!ComputeResult(two, '*', xn, temp3)) {
            printf("error 4\n");
            return;
        }

        if (!ComputeResult(temp1, '/', temp3, xn1)) {
            printf("error 5\n");
            return;
        }

        TermToAscii(xn1, buffer, NORMAL);
        printf("%2d %s\n", i, buffer);

        /* Are we done ? */
        if (!NormAbsCmp(xn, xn1))
            break;

        if (i > 49) {
            printf("\nLast entry is the closest approxima"
                   "tion. Algorithm no longer converging\n");
            break;
        }

        /* The current xn1 becomes the next xn */
        TermCopy(xn, xn1);
    }

    /* print the square root */
    TermToAscii(xn1, buffer, NORMAL);
    printf("\nSquare root: %s\n", buffer);

    /* print the square of the square root */
    TermCopy(xn, xn1);
    if (!ComputeResult(xn, '*', xn1, product)) {
        printf("error 6\n");
        return;
    }

    TermToAscii(product, buffer, NORMAL);
    printf("Computed square: %s\n", buffer);

    /* print the difference */
    if (!ComputeResult(square, '-', product, xn)) {
        printf("error 7\n");
        return;
    }
    TermToAscii(xn, buffer, NORMAL);

    printf("Delta: %s\n", buffer);

    /* now test the built-in floating-point routines */
    fsquare = atof(argv[1]);
    fsqrt = sqrt(fsquare);

    printf("\n\nMath lib root: %4.12f, square: %4.12f\n",
           fsqrt, fsqrt * fsqrt);

    return;
}
